import math
import numpy as np
from scipy import signal
from PIL import Image, ImageDraw
from scipy.ndimage import gaussian_filter
import os

Part 1

Question 1-2

import ncc

def MakeGaussianPyramind(image, scale, minsize):
    '''
    purpose: creates a Gaussian pyramid for an image
    inputs: image - jpg, scale - factor by which to reduce image size, minsize - min dimension of reduced image
    outputs: list of [original_image, PIL images of reduced size]
    '''
    # generate appropriate sigma
    sigma=1/(2*scale)
    # open image
    im = Image.open(image)
    pyr = []
    pyr.append(np.float32(im))
    # set while loop with min dim as stop condition
    x, y = im.width, im.height
    min_dim = np.max([x, y])
    while int(min_dim*0.75) > minsize: # check condition
        if len(np.float32(im).shape) < 3: # if grayscale, do not split channels for separate filtering
            im2 = Image.fromarray(gaussian_filter(im, sigma))
            # resize image
            imnew = im2.resize((int(x*scale), int(y*scale)), Image.BICUBIC)
            # convert back to np array
            resized = np.float32(imnew)
            # set new ref image as recently downsampled one
            im = imnew
        else: # if RGB. apply filters separately
            r, g, b = im.split() # if multi-channel, split channels and filter before reassmbeling
            r, g, b = Image.fromarray(gaussian_filter(r, sigma)), Image.fromarray(gaussian_filter(g, sigma)), Image.fromarray(gaussian_filter(b, sigma))
            # clip values and merge channels
            im2 = Image.merge('RGB', (r, g, b))
            # resize image
            imnew = im2.resize((int(x*scale), int(y*scale)), Image.BICUBIC)
            # convert back to numpy arrays 
            r, g, b = imnew.split()
            tmp = Image.merge('RGB', (r, g, b))
            resized = np.float32(tmp)
            # set new ref image as recently downsampled one
            im = imnew
        pyr.append(resized)
        # update min_dim
        x, y = imnew.width, imnew.height
        min_dim = np.max([x, y])
    return pyr
        

print(os.getcwd())
filepath=r"hw2part2\faces\judybats.jpg"
pyr = MakeGaussianPyramind(image=filepath, scale=0.75, minsize=50)
print('done')

c:\Users\Sadie\Documents\GitHub\CPSC_425\A2
done

Question 3

def ShowGaussianPyramid(pyr):
    '''
    purpose: display Gaussian pyramid
    inputs: list of np arrays
    outputs: None (display all images in pyramid side by side)
    '''
    heights = list(map(lambda x: x.shape[0], pyr)) # get dims of all imgs in list to know dims needed for template image
    widths = list(map(lambda x: x.shape[1], pyr))
    w_tot, h_tot = np.sum(widths), np.max(heights)
    image = Image.new(mode="L", size=(w_tot, h_tot), color="#ffff")
    offset_x, offset_y = 0, 0
    for im in pyr:
        im_pil = Image.fromarray(im)
        image.paste(im_pil,(offset_x, offset_y))
        offset_x += im_pil.width
        #offset_y = im_pil.height
    print(image.size)
    display(image)

pyr = MakeGaussianPyramind(image=filepath, scale=0.75, minsize=100)
ShowGaussianPyramid(pyr)

(1089, 342)

Question 4

def FindTemplate(pyramid, template, threshold):
    '''
    purpose: find and mark regions in set of images where template matches above a certain threshold
    inputs: list of images (pyramid), template, threshold for matching
    outputs: annotated image
    '''
    # reduce template size
    temp_w = 15
    ar = template.height / template.width # keep same aspect ratio
    template = template.resize((temp_w,int(temp_w*ar)), Image.BICUBIC)
    # loop through images
    first_im = Image.fromarray(pyramid[0]).convert('RGB')
    colors=['red','orange','yellow','green','blue','purple','pink','black']
    for i, im in enumerate(pyramid):
        res = ncc.normxcorr2D(im, template)
        y, x = np.where(res > threshold)
        for j in range(len(x)):
            def rescale_x(x):
                ''' 
                purpose: rescale a x coordinate from a smaller image to correspond to the same location in a larger image
                input: x, the coordinate
                '''
                return (x/Image.fromarray(im).width)*first_im.width
            
            def rescale_y(y):
                ''' 
                purpose: rescale a y coordinate from a smaller image to correspond to the same location in a larger image
                input: y, the coordinate
                '''
                return (y/Image.fromarray(im).height)*first_im.height
            # create corners of box to be drawn
            temp_w, temp_h = template.width, template.height
            xc,yc = rescale_x(x[j]),rescale_y(y[j]) #,rescale_x(x2),rescale_y(y2)
            x1, x2 = xc - (temp_w/0.75**i)//2, xc + (temp_w/0.75**i)//2
            y1, y2 = yc - (temp_h/0.75**i)//2, yc + (temp_h/0.75**i)//2
            # rescale coordinates of high correlation for plotting on original image
            
            draw = ImageDraw.Draw(first_im)
            draw.line((x1,y1,x2,y1),fill=colors[i],width=2)
            draw.line((x1,y1,x1,y2),fill=colors[i],width=2)
            draw.line((x1,y2,x2,y2),fill=colors[i],width=2)
            draw.line((x2,y1,x2,y2),fill=colors[i],width=2)
            #del draw
    #display(first_im)
    #first_im.show()
    first_im.save("template_match.tif")
    return first_im
            

temp = Image.open(r"hw2part2\faces\template.jpg")
FindTemplate(pyr, temp, 0.72)

c:\Users\Sadie\Documents\GitHub\CPSC_425\A2\ncc.py:59: RuntimeWarning: divide by zero encountered in divide
  nxcorr = np.where(denom < tol, 0, numer/denom)

Question 5

# testing thresholds
thresh=0.8

filepath=r"hw2part2\faces\judybats.jpg"
pyr = MakeGaussianPyramind(image=filepath, scale=0.75, minsize=50)
x = FindTemplate(pyr, temp, thresh)
filepath=r"hw2part2\faces\students.jpg"
pyr = MakeGaussianPyramind(image=filepath, scale=0.75, minsize=50)
x = FindTemplate(pyr, temp, thresh)
filepath=r"hw2part2\faces\sports.jpg"
pyr = MakeGaussianPyramind(image=filepath, scale=0.75, minsize=50)
x = FindTemplate(pyr, temp, thresh)
filepath=r"hw2part2\faces\family.jpg"
pyr = MakeGaussianPyramind(image=filepath, scale=0.75, minsize=50)
x = FindTemplate(pyr, temp, thresh)

thresh=0.4

filepath=r"hw2part2\faces\tree.jpg"
pyr = MakeGaussianPyramind(image=filepath, scale=0.75, minsize=50)
x = FindTemplate(pyr, temp, thresh)
filepath=r"hw2part2\faces\students.jpg"
pyr = MakeGaussianPyramind(image=filepath, scale=0.75, minsize=50)
x = FindTemplate(pyr, temp, thresh)
filepath=r"hw2part2\faces\sports.jpg"
pyr = MakeGaussianPyramind(image=filepath, scale=0.75, minsize=50)
x = FindTemplate(pyr, temp, thresh)
filepath=r"hw2part2\faces\family.jpg"
pyr = MakeGaussianPyramind(image=filepath, scale=0.75, minsize=50)
x = FindTemplate(pyr, temp, thresh)

thresh=0.65

filepath=r"hw2part2\faces\judybats.jpg"
pyr = MakeGaussianPyramind(image=filepath, scale=0.75, minsize=50)
x = FindTemplate(pyr, temp, thresh)
filepath=r"hw2part2\faces\students.jpg"
pyr = MakeGaussianPyramind(image=filepath, scale=0.75, minsize=50)
x = FindTemplate(pyr, temp, thresh)
filepath=r"hw2part2\faces\tree.jpg"
pyr = MakeGaussianPyramind(image=filepath, scale=0.75, minsize=50)
x = FindTemplate(pyr, temp, thresh)
filepath=r"hw2part2\faces\fans.jpg"
pyr = MakeGaussianPyramind(image=filepath, scale=0.75, minsize=50)
x = FindTemplate(pyr, temp, thresh)

# best threshold
thresh=0.58

# judybats
filepath=r"hw2part2\faces\judybats.jpg"
pyr = MakeGaussianPyramind(image=filepath, scale=0.75, minsize=50)
x = FindTemplate(pyr, temp, thresh)
display(x)

# False positives: 1
# False negatives: 1

# students
filepath=r"hw2part2\faces\students.jpg"
pyr = MakeGaussianPyramind(image=filepath, scale=0.75, minsize=50)
#ShowGaussianPyramid(pyr)
x = FindTemplate(pyr, temp, thresh)
display(x)

# False positives: 3
# False negatives: 4

filepath=r"hw2part2\faces\tree.jpg"
pyr = MakeGaussianPyramind(image=filepath, scale=0.75, minsize=50)
x = FindTemplate(pyr, temp, thresh)
display(x)
# False positives: 2
# False negatives: 0

filepath=r"hw2part2\faces\family.jpg"
pyr = MakeGaussianPyramind(image=filepath, scale=0.75, minsize=50)
#ShowGaussianPyramid(pyr)
x = FindTemplate(pyr, temp, thresh)
display(x)
# False positives: 0
# False negatives: 1

filepath=r"hw2part2\faces\fans.jpg"
pyr = MakeGaussianPyramind(image=filepath, scale=0.75, minsize=50)
x = FindTemplate(pyr, temp, thresh)
display(x)
# False positives: 3
# False negatives: 3

filepath=r"hw2part2\faces\sports.jpg"
pyr = MakeGaussianPyramind(image=filepath, scale=0.75, minsize=50)
x = FindTemplate(pyr, temp, thresh)
display(x)
# False positives: 1
# False negatives: 1

A threshold of 0.58 achieves an approximately equal error rate across all tested images. While there is an approximately equal error rate across all images, this threshold shows the best performance for the 'students' image, where nearly all faces are correctly detected with only a small fraction of false negatives. Higher thresholds led to missed faces in several of the images, whereas lower thresholds led to a number of false positives, such as in the image of the tree where no true faces are present.

Question 6

Using a threshold of 0.58 on the NCC, the recall rates on each image are as follows:

import pandas as pd 

tps = [4,23,0,2,0,0]
pos = [5,27,0,3,3,1]
recall = []
for i in range(len(pos)):
    if pos[i] != 0:
        recall.append(tps[i]/pos[i])
    else:
        recall.append('undefined')

summary = {"Image": ['judybats', 'students', 'tree', 'family', 'fans', 'sports'], 
          "True Positives":tps, 
          "Positives":pos, 
          "Recall": recall
}
print(pd.DataFrame(summary))

      Image  True Positives  Positives     Recall
0  judybats               4          5        0.8
1  students              23         27   0.851852
2      tree               0          0  undefined
3    family               2          3   0.666667
4      fans               0          3        0.0
5    sports               0          1        0.0

The recall rate using NCC is rather low in particular for images that have very few actual positives for faces. For example, in the 'family' photo, while only 1 face was missed, the recall was only 66%. Thus recall seems like it is only a useful metric if there is a sizeable number of potential detections. Recall is also problematic as it does not report at all on false positives, which could be abundant amidst the correctly identified true positives, so this should be kept in mind as well.

Part 2: Image Blending

def MakeLaplacianPyramind(image, scale, minsize):
    '''
    purpose: creates a Laplacian pyramid for an image
    inputs: image - jpg, scale - factor by which to reduce image size, minsize - min dimension of reduced image
    outputs: list of [original_image, PIL images of reduced size]
    '''
    # initialize sigma
    sigma=1/(2*scale)
    # generate Gaussian pyramid
    Gauss_pyr = MakeGaussianPyramind(image, scale, minsize)
    Laplace_pyr = []
    for i, img in enumerate(Gauss_pyr):
        r, g, b = img[:,:,0], img[:,:,1], img[:,:,2]
        r2, g2, b2 = gaussian_filter(r, sigma), gaussian_filter(g, sigma), gaussian_filter(b, sigma)
        # clip to 0-255 and convert to image
        def clip_to(ch, min=0, max=255):
            return Image.fromarray(np.clip(ch, min, max).astype(np.uint8))
        r2_pil, g2_pil, b2_pil = clip_to(r2), clip_to(g2), clip_to(b2)
        # create smoothed image
        smoothed = Image.merge('RGB', (r2_pil, g2_pil, b2_pil))
        # if at end of pyramid, append tiny Gaussian
        if i == len(Gauss_pyr)-1:
            Laplace_pyr.append(np.float32(smoothed))
        # otherwise compute Laplacian as difference between original and blurred
        else:
            # subtract each channel
            r_dif, g_dif, b_dif = r-r2, g-g2, b-b2
            Laplace_img =  np.stack(np.array([r_dif, g_dif, b_dif]), axis=2)
            Laplace_pyr.append(np.float32(Laplace_img))
    return Laplace_pyr

Question 3

def ShowLaplacianPyramid(pyr):
    '''
    purpose: display Laplacian pyramid
    inputs: list of np arrays
    outputs: None (display all images in pyramid side by side)
    '''
    heights = list(map(lambda x: x.shape[0], pyr)) # get dims of all imgs in list to know dims needed for template image
    widths = list(map(lambda x: x.shape[1], pyr))
    w_tot, h_tot = np.sum(widths), np.max(heights)
    image = Image.new(mode="RGB", size=(w_tot, h_tot), color="#ffff")
    offset_x, offset_y = 0, 0
    for im in pyr:
        # adjust pixel values for viewing and discretize
        im += 128
        # normalize pixel range to 0-255 for improved contrast
        im = (im - np.min(im))/(np.max(im)-np.min(im))*255
        # clip
        im = np.clip(im, 0, 255).astype(np.uint8)
        # convert to RGB PIL Image
        r,g,b = Image.fromarray(im[:,:,0]), Image.fromarray(im[:,:,1]), Image.fromarray(im[:,:,2])
        im_pil = Image.merge('RGB', (r,g,b))
        image.paste(im_pil,(offset_x, offset_y))
        offset_x += im_pil.width
    display(image)

def ShowGaussianPyramidRGB(pyr):
    '''
    purpose: display Laplacian pyramid
    inputs: list of np arrays
    outputs: None (display all images in pyramid side by side)
    '''
    heights = list(map(lambda x: x.shape[0], pyr)) # get dims of all imgs in list to know dims needed for template image
    widths = list(map(lambda x: x.shape[1], pyr))
    w_tot, h_tot = np.sum(widths), np.max(heights)
    image = Image.new(mode="RGB", size=(w_tot, h_tot), color="#ffff")
    offset_x, offset_y = 0, 0
    for im in pyr:
        # adjust pixel values for viewing and discretize
        im = im.astype(np.uint8)
        # convert to RGB PIL Image
        r,g,b = Image.fromarray(im[:,:,0]), Image.fromarray(im[:,:,1]), Image.fromarray(im[:,:,2])
        im_pil = Image.merge('RGB', (r,g,b))
        image.paste(im_pil,(offset_x, offset_y))
        offset_x += im_pil.width
    display(image)

scale=0.75
minsize=50

# violet
filepath=r"hw2part2\violet.jpg"
violet_g = MakeGaussianPyramind(image=filepath, scale=scale, minsize=minsize)
ShowGaussianPyramidRGB(violet_g)
violet = MakeLaplacianPyramind(image=filepath, scale=scale, minsize=minsize)
ShowLaplacianPyramid(violet)

# orchid
filepath=r"hw2part2\orchid.jpg"
orchid_g = MakeGaussianPyramind(image=filepath, scale=scale, minsize=minsize)
ShowGaussianPyramidRGB(orchid_g)
orchid = MakeLaplacianPyramind(image=filepath, scale=scale, minsize=minsize)
ShowLaplacianPyramid(orchid)

Question 4

def ReconstructGaussianFromLaplacianPyramid(lPyramid):
    '''
    purpose: takes a Laplacian pyramid and reconstructs the Gaussian pyramid of that image from a Laplacian on
    input: lPyramid, a Laplacian pyramid
    output: gPyramid, a Gaussian pyramid
    '''
    gPyramid = []
    # add first entry to gPyramid since its the same as lPyramid
    gPyramid.append(lPyramid[-1])
    for i in range(1,len(lPyramid)):
        curr = gPyramid[-1]
        next = lPyramid[-i-1]
        # get size of next element in the pyr
        next_w, next_h = int(next.shape[1]), int(next.shape[0])
        # convert Gaussian to Image format for resizing
        r, g, b = Image.fromarray(curr[:,:,0].astype(np.uint8)), Image.fromarray(curr[:,:,1].astype(np.uint8)), Image.fromarray(curr[:,:,2].astype(np.uint8))
        curr_PIL = Image.merge('RGB', (r, g, b))
        # upsample/resize the image
        up_PIL = curr_PIL.resize((next_w, next_h), Image.BICUBIC)
        up = np.float32(up_PIL)
        # add details back using laplacian
        g_img = up + next
        # append reconstruction to gPyramid
        gPyramid.append(g_img)
    # reverse order of pyramid
    gPyramid.reverse()
    return gPyramid

Small Artifacts

# violet
filepath=r"hw2part2\violet.jpg"
violet_g1 = MakeGaussianPyramind(image=filepath, scale=scale, minsize=minsize)
violet_l = MakeLaplacianPyramind(image=filepath, scale=scale, minsize=minsize)
violet_g2 = ReconstructGaussianFromLaplacianPyramid(violet_l)
print('original')
ShowGaussianPyramidRGB(violet_g1)
print('reconstructed')
ShowGaussianPyramidRGB(violet_g2)

# orchid
filepath=r"hw2part2\orchid.jpg"
orchid_g1 = MakeGaussianPyramind(image=filepath, scale=scale, minsize=minsize)
orchid_l = MakeLaplacianPyramind(image=filepath, scale=scale, minsize=minsize)
orchid_g2 = ReconstructGaussianFromLaplacianPyramid(orchid_l)
print('original')
ShowGaussianPyramidRGB(orchid_g1)
print('reconstructed')
ShowGaussianPyramidRGB(orchid_g2)

original

reconstructed

original

reconstructed

Artifacts in the merged image due to aliasing can be mitigated by decreasing the sigma used for the blurring, however, this also reduces the integrity of the resulting reproduction (e.g., edges are less sharp, colours are less saturated).

# edit Gaussian pyramid function to reduce artifacts
# specifically, increase sigma from 1/2*scale to 1/2.3*scale
def MakeGaussianPyramind_noartifacts(image, scale, minsize):
    '''
    purpose: creates a Gaussian pyramid for an image
    inputs: image - jpg, scale - factor by which to reduce image size, minsize - min dimension of reduced image
    outputs: list of [original_image, PIL images of reduced size]
    '''
    # generate appropriate sigma
    sigma=1/(2.3*scale)
    # open image
    im = Image.open(image)
    pyr = []
    pyr.append(np.float32(im))
    # set while loop with min dim as stop condition
    x, y = im.width, im.height
    min_dim = np.max([x, y])
    while int(min_dim*0.75) > minsize: # check condition
        if len(np.float32(im).shape) < 3: # if grayscale, do not split channels for separate filtering
            im2 = Image.fromarray(gaussian_filter(im, sigma))
            # resize image
            imnew = im2.resize((int(x*scale), int(y*scale)), Image.BICUBIC)
            # convert back to np array
            resized = np.float32(imnew)
            # set new ref image as recently downsampled one
            im = imnew
        else: # if RGB. apply filters separately
            r, g, b = im.split() # if multi-channel, split channels and filter before reassmbeling
            r, g, b = Image.fromarray(gaussian_filter(r, sigma)), Image.fromarray(gaussian_filter(g, sigma)), Image.fromarray(gaussian_filter(b, sigma))
            # clip values and merge channels
            im2 = Image.merge('RGB', (r, g, b))
            # resize image
            imnew = im2.resize((int(x*scale), int(y*scale)), Image.BICUBIC)
            # convert back to numpy arrays 
            r, g, b = imnew.split()
            tmp = Image.merge('RGB', (r, g, b))
            resized = np.float32(tmp)
            # set new ref image as recently downsampled one
            im = imnew
        pyr.append(resized)
        # update min_dim
        x, y = imnew.width, imnew.height
        min_dim = np.max([x, y])
    return pyr

def MakeLaplacianPyramind_noartifacts(image, scale, minsize):
    '''
    purpose: creates a Laplacian pyramid for an image
    inputs: image - jpg, scale - factor by which to reduce image size, minsize - min dimension of reduced image
    outputs: list of [original_image, PIL images of reduced size]
    '''
    # initialize sigma
    sigma=1/(2.3*scale)
    # generate Gaussian pyramid
    Gauss_pyr = MakeGaussianPyramind(image, scale, minsize)
    Laplace_pyr = []
    for i, img in enumerate(Gauss_pyr):
        r, g, b = img[:,:,0], img[:,:,1], img[:,:,2]
        r2, g2, b2 = gaussian_filter(r, sigma), gaussian_filter(g, sigma), gaussian_filter(b, sigma)
        # clip to 0-255 and convert to image
        def clip_to(ch, min=0, max=255):
            return Image.fromarray(np.clip(ch, min, max).astype(np.uint8))
        r2_pil, g2_pil, b2_pil = clip_to(r2), clip_to(g2), clip_to(b2)
        # create smoothed image
        smoothed = Image.merge('RGB', (r2_pil, g2_pil, b2_pil))
        # if at end of pyramid, append tiny Gaussian
        if i == len(Gauss_pyr)-1:
            Laplace_pyr.append(np.float32(smoothed))
        # otherwise compute Laplacian as difference between original and blurred
        else:
            # subtract each channel
            r_dif, g_dif, b_dif = r-r2, g-g2, b-b2
            Laplace_img =  np.stack(np.array([r_dif, g_dif, b_dif]), axis=2)
            Laplace_pyr.append(np.float32(Laplace_img))
    return Laplace_pyr

# violet
filepath=r"hw2part2\violet.jpg"
violet_l = MakeLaplacianPyramind_noartifacts(image=filepath, scale=scale, minsize=minsize)
violet_g2 = ReconstructGaussianFromLaplacianPyramid(violet_l)
print('original')
ShowGaussianPyramidRGB(violet_g1)
print('reconstructed-no artifacts')
ShowGaussianPyramidRGB(violet_g2)

# orchid
filepath=r"hw2part2\orchid.jpg"
orchid_l = MakeLaplacianPyramind_noartifacts(image=filepath, scale=scale, minsize=minsize)
orchid_g2 = ReconstructGaussianFromLaplacianPyramid(orchid_l)
print('original')
ShowGaussianPyramidRGB(orchid_g1)
print('reconstructed-less artifacts')
ShowGaussianPyramidRGB(orchid_g2)

original

reconstructed-no artifacts

original

reconstructed-less artifacts

Question 5

filepath=r"hw2part2\orchid_mask.bmp"
gaussianM = MakeGaussianPyramind(image=filepath, scale=scale, minsize=minsize)
ShowGaussianPyramid(gaussianM)

(1840, 341)

Question 6

Despite the adjustment of the blurring prior to reconstruction, there is still some small artifacting in the resulting reconstructions. This seems to be a necessary trade off in order to preserve the edges and saturation of the original component images.

# Compose a single Laplacian pyramid from two, level by level
lapA = orchid_l
lapB = violet_l
compLaplacian = []
for i, lvl in enumerate(gaussianM):
    # reshape binary mask for compatibility with RGB
    mask = gaussianM[i].reshape(gaussianM[i].shape[0], gaussianM[i].shape[1], 1)
    entry = lapA[i] * mask/255 +  lapB[i] * (1 - mask/255)
    compLaplacian.append(entry)
# Reconstruct Gaussian from the composed Laplacian
recon=ReconstructGaussianFromLaplacianPyramid(compLaplacian)

# display the highest resolution image
curr = recon[0]
r, g, b = Image.fromarray(curr[:,:,0].astype(np.uint8)), Image.fromarray(curr[:,:,1].astype(np.uint8)), Image.fromarray(curr[:,:,2].astype(np.uint8))
res = Image.merge('RGB', (r, g, b))
display(res)

# once again, adjust sigma to reduce artifacts

def MakeGaussianPyramind_noartifacts(image, scale, minsize):
    '''
    purpose: creates a Gaussian pyramid for an image
    inputs: image - jpg, scale - factor by which to reduce image size, minsize - min dimension of reduced image
    outputs: list of [original_image, PIL images of reduced size]
    '''
    # generate appropriate sigma
    sigma=1/(2.6*scale)
    # open image
    im = Image.open(image)
    pyr = []
    pyr.append(np.float32(im))
    # set while loop with min dim as stop condition
    x, y = im.width, im.height
    min_dim = np.max([x, y])
    while int(min_dim*0.75) > minsize: # check condition
        if len(np.float32(im).shape) < 3: # if grayscale, do not split channels for separate filtering
            im2 = Image.fromarray(gaussian_filter(im, sigma))
            # resize image
            imnew = im2.resize((int(x*scale), int(y*scale)), Image.BICUBIC)
            # convert back to np array
            resized = np.float32(imnew)
            # set new ref image as recently downsampled one
            im = imnew
        else: # if RGB. apply filters separately
            r, g, b = im.split() # if multi-channel, split channels and filter before reassmbeling
            r, g, b = Image.fromarray(gaussian_filter(r, sigma)), Image.fromarray(gaussian_filter(g, sigma)), Image.fromarray(gaussian_filter(b, sigma))
            # clip values and merge channels
            im2 = Image.merge('RGB', (r, g, b))
            # resize image
            imnew = im2.resize((int(x*scale), int(y*scale)), Image.BICUBIC)
            # convert back to numpy arrays 
            r, g, b = imnew.split()
            tmp = Image.merge('RGB', (r, g, b))
            resized = np.float32(tmp)
            # set new ref image as recently downsampled one
            im = imnew
        pyr.append(resized)
        # update min_dim
        x, y = imnew.width, imnew.height
        min_dim = np.max([x, y])
    return pyr

def MakeLaplacianPyramind_noartifacts(image, scale, minsize):
    '''
    purpose: creates a Laplacian pyramid for an image
    inputs: image - jpg, scale - factor by which to reduce image size, minsize - min dimension of reduced image
    outputs: list of [original_image, PIL images of reduced size]
    '''
    # initialize sigma
    sigma=1/(2.6*scale)
    # generate Gaussian pyramid
    Gauss_pyr = MakeGaussianPyramind(image, scale, minsize)
    Laplace_pyr = []
    for i, img in enumerate(Gauss_pyr):
        r, g, b = img[:,:,0], img[:,:,1], img[:,:,2]
        r2, g2, b2 = gaussian_filter(r, sigma), gaussian_filter(g, sigma), gaussian_filter(b, sigma)
        # clip to 0-255 and convert to image
        def clip_to(ch, min=0, max=255):
            return Image.fromarray(np.clip(ch, min, max).astype(np.uint8))
        r2_pil, g2_pil, b2_pil = clip_to(r2), clip_to(g2), clip_to(b2)
        # create smoothed image
        smoothed = Image.merge('RGB', (r2_pil, g2_pil, b2_pil))
        # if at end of pyramid, append tiny Gaussian
        if i == len(Gauss_pyr)-1:
            Laplace_pyr.append(np.float32(smoothed))
        # otherwise compute Laplacian as difference between original and blurred
        else:
            # subtract each channel
            r_dif, g_dif, b_dif = r-r2, g-g2, b-b2
            Laplace_img =  np.stack(np.array([r_dif, g_dif, b_dif]), axis=2)
            Laplace_pyr.append(np.float32(Laplace_img))
    return Laplace_pyr

Cup

filepath=r"hw2part2\blue_cup.jpg"
lapB = MakeLaplacianPyramind_noartifacts(image=filepath, scale=scale, minsize=minsize)
filepath=r"hw2part2\green_cup.jpg"
lapA = MakeLaplacianPyramind_noartifacts(image=filepath, scale=scale, minsize=minsize)
# create mask
filepath=r"hw2part2\cup_mask.bmp"
gaussianM = MakeGaussianPyramind_noartifacts(image=filepath, scale=scale, minsize=minsize)

# Compose a single Laplacian pyramid from two, level by level
compLaplacian = []
for i, lvl in enumerate(gaussianM):
    # reshape binary mask for compatibility with RGB
    mask = gaussianM[i].reshape(gaussianM[i].shape[0], gaussianM[i].shape[1], 1)
    entry = lapA[i] * mask/255 +  lapB[i] * (1 - mask/255)
    compLaplacian.append(entry)
# Reconstruct Gaussian from the composed Laplacian
recon=ReconstructGaussianFromLaplacianPyramid(compLaplacian)
#ShowGaussianPyramidRGB(recon)

# display the highest resolution image
curr = recon[0]
r, g, b = Image.fromarray(curr[:,:,0].astype(np.uint8)), Image.fromarray(curr[:,:,1].astype(np.uint8)), Image.fromarray(curr[:,:,2].astype(np.uint8))
res = Image.merge('RGB', (r, g, b))
display(res)

Fruit

filepath=r"hw2part2\tomato.jpg"
lapA = MakeLaplacianPyramind_noartifacts(image=filepath, scale=scale, minsize=minsize)
filepath=r"hw2part2\apple.jpg"
lapB = MakeLaplacianPyramind_noartifacts(image=filepath, scale=scale, minsize=minsize)
# create mask
filepath=r"hw2part2\tomato_mask.bmp"
gaussianM = MakeGaussianPyramind_noartifacts(image=filepath, scale=scale, minsize=minsize)

# Compose a single Laplacian pyramid from two, level by level
compLaplacian = []
for i, lvl in enumerate(gaussianM):
    # reshape binary mask for compatibility with RGB
    mask = gaussianM[i].reshape(gaussianM[i].shape[0], gaussianM[i].shape[1], 1)
    entry = lapA[i] * mask/255 +  lapB[i] * (1 - mask/255)
    compLaplacian.append(entry)
# Reconstruct Gaussian from the composed Laplacian
recon=ReconstructGaussianFromLaplacianPyramid(compLaplacian)
#ShowGaussianPyramidRGB(recon)

# display the highest resolution image
curr = recon[0]
r, g, b = Image.fromarray(curr[:,:,0].astype(np.uint8)), Image.fromarray(curr[:,:,1].astype(np.uint8)), Image.fromarray(curr[:,:,2].astype(np.uint8))
res = Image.merge('RGB', (r, g, b))
display(res)